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We investigate the eigenvalue spectrum of the staggered Dirac matrix 
in SU(3) gauge theory and in full QCD as well as in quenched U(l) theory 
on various lattice sizes. As a measure of the fluctuation properties of the 
eigenvalues, we study the nearest-neighbor spacing distribution, P(s). We 
further study lattice QCD at nonzero chemical potential, fi, by construct- 
ing the spacing distribution of adjacent eigenvalues in the complex plane. 
We find that in all regions of their phase diagrams, compact lattice gauge 
theories have bulk spectral correlations given by random matrix theory, 
which is an indication for quantum chaos. 

1. QCD at Nonzero Temperature 

The properties of the eigenvalues of the Dirac operator are of great im- 
portance for the understanding of certain features of QCD. For example, 
the accumulation of small eigenvalues is, via the Banks-Casher formula [1], 
related to the spontaneous breaking of chiral symmetry. Recently, the fluc- 
tuation properties of the eigenvalues in the bulk of the spectrum have also 
attracted attention. It was shown in Ref. [2] that on the scale of the mean 
level spacing they are described by random matrix theory (RMT). In partic- 
ular, the nearest-neighbor spacing distribution P(s), i.e., the distribution of 
spacings s between adjacent eigenvalues on the unfolded scale, agrees with 
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the Wigner surmise of RMT. According to the Bohigas-Giannoni-Schmit 
conjecture [3], quantum systems whose classical counterparts are chaotic 
have a nearest-neighbor spacing distribution given by RMT whereas sys- 
tems whose classical counterparts are integrable obey a Poisson distribu- 
tion, Pp(s) = e~ s . Therefore, the specific form of P(s) is often taken as a 
criterion for the presence or absence of "quantum chaos" . 

In RMT, one has to distinguish several universality classes which are 
determined by the symmetries of the system. For the case of the QCD 
Dirac operator, this classification was done in Ref. [4]. Depending on the 
number of colors and the representation of the quarks, the Dirac opera- 
tor is described by one of the three chiral ensembles of RMT. As far as 
the fluctuation properties in the bulk of the spectrum are concerned, the 
predictions of the chiral ensembles are identical to those of the ordinary 
ensembles [5]. In Ref. [2], the Dirac matrix was studied in SU(2) using both 
staggered and Wilson fermions which correspond to the chiral symplec- 
tic and orthogonal ensemble, respectively. Here [6], we study SU(3) with 
staggered fermions which corresponds to the chiral unitary ensemble. The 
RMT result for the nearest-neighbor spacing distribution can be expressed 
in terms of so-called prolate spheroidal functions, see Ref. [7]. A very good 
approximation to P(s) is provided by the Wigner surmise for the unitary 
ensemble, 

JVM = - 2 s\-^/* . (i) 

We generated gauge field configurations using the standard Wilson pla- 
quette action for SU(3) with and without dynamical fermions in the Kogut- 
Susskind prescription. We have worked on a 6 3 x 4 lattice with various val- 
ues of the inverse gauge coupling j3 = Q/g 2 both in the confinement and 
deconfinement phase. We typically produced 10 independent equilibrium 
configurations for each (3. Because of the spectral ergodicity property of 
RMT one can replace ensemble averages by spectral averages if one is only 
interested in bulk properties. 

The Dirac operator,^) = is anti-hermitian so that all eigenvalues 

are imaginary. For convenience, we denote them by i\n and refer to trie X n 
as the eigenvalues. Because of {^,75} = the nonzero occur in pairs 
of opposite sign. All spectra were checked against the analytical sum rules 
J2n A ra = and J2x n >o = 3V^, where V is the lattice volume. To construct 
the nearest-neighbor spacing distribution from the eigenvalues, one first has 
to "unfold" the spectra [7]. 

Figure 1 compares P(s) of full QCD with Nf = 3 flavors and quark 
mass ma = 0.05 to the RMT result. In the confinement as well as in the 
deconfinement phase we observe agreement with RMT up to very high (3 
(not shown). The observation that P(s) is not influenced by the presence 
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Figure 1. The nearest-neighbor spacing distribution P(s) on a 6 3 x 4 lattice in full QCD 
(histograms) compared with the random matrix result (solid lines) . There are no changes 
in P(s) across the deconfinement phase transition. 



of dynamicaf quarks could have been expected from the results of Ref. [5], 
which apply to the case of massless quarks. Our results, and those of [2], 
indicate that massive dynamical quarks do not affect P(s) either. 

No signs for a transition to Poisson regularity are found. The decon- 
finement phase transition does not seem to coincide with a transition in 
the spacing distribution. For very large values of (3 far into the deconfine- 
ment region, the eigenvalues start to approach the degenerate eigenvalues 
of the free theory, given by A 2 = X^=i s in 2 (27rn At /L M )/a 2 , where a is the 
lattice constant, L M is the number of lattice sites in the //-direction, and 
n M = 0, . . . , Lfj, — 1. In this case, the spacing distribution is neither Wigner 
nor Poisson. It is possible to lift the degeneracies of the free eigenvalues 
using an asymmetric lattice where L x , L y , etc. are relative primes and, for 
large lattices, the distribution is then Poisson, Pp(s) = e~ s , see Fig. 2. 

2. QCD at Nonzero Density 

Physical systems which are described by non-hermitian operators have at- 
tracted a lot of attention recently, among others QCD at nonzero chemical 
potential fi [8]. There, the Dirac operator loses its hermiticity properties 
so that its eigenvalues become complex. The aim of the present analysis is 
to investigate whether non-hermitian RMT is able to describe the fluctua- 
tion properties of the complex eigenvalues of the QCD Dirac operator. The 
eigenvalues are generated on the lattice for various values of \x. We apply 
a two-dimensional unfolding procedure to separate the average eigenvalue 
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Figure 2. Nearest-neighbor spacing distribution P(s) for the free Dirac operator on a 
53 x 47 x 43 x 41 lattice compared with a Poisson distribution, e _s . 



density from the fluctuations and construct the nearest-neighbor spacing 
distribution, P(s), of adjacent eigenvalues in the complex plane. The data 
are then compared to analytical predictions of non-hermitian RMT. 

We start with a few definitions. A formulation of the QCD Dirac oper- 
ator at \i / on the lattice in the staggered scheme is given by [9] 

M x>y (U,n) = — [Uu{x)r] u {x)5y^ - h.c] 

v=x,y,z 

+ Ya [^K^^W " UHVMV)*-"^] ( 2 ) 

with the link variables U and the staggered phases rj. 

We consider the gauge group SU(3) which corresponds to the symmetry 
class of the chiral unitary ensemble of RMT [4, 10]. At zero chemical poten- 
tial, all Dirac eigenvalues are purely imaginary, and the nearest-neighbor 
spacing distribution, P(s), of the lattice data agrees with the Wigner sur- 
mise of hermitian RMT, Eq. (1), both in the confinement and in the de- 
confinement phase (see Fig. 1). This finding implies strong correlations of 
the eigenvalues. For uncorrelated eigenvalues P(s) is given by the Poisson 
distribution. 

For a complex spectrum, we define P(s) to represent the spacing distri- 
bution of nearest neighbors in the complex plane, i.e., for each eigenvalue 
zq one identifies the eigenvalue z\ for which s = \zi — zo\ is a minimum 
[11]. After ensemble averaging, one obtains a function P(s,zq) which, in 
general, depends on zq. The dependence on zq can be eliminated by un- 
folding the spectrum, i.e., by applying a local rescaling of the energy scale 
so that the average spectral density is constant in a bounded region in the 
complex plane and zero outside [12]. After unfolding, a spectral average 
over zq yields P(s). 
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For n > 0, the eigenvalues of the matrix in Eq. (2) move into the 
complex plane. If the real and imaginary parts of the strongly correlated 
eigenvalues have approximately the same average magnitude, the system 
should be described by the Ginibre ensemble of non-hermitian RMT [13]. 
In the Ginibre ensemble, the average spectral density is already constant 
inside a circle and zero outside, respectively. In this case, unfolding is not 
necessary, and P(s) is given by [11] 



where e n (x) = Em=o x "7 m! and c = J™ dss P( s ) = 1-1429.... This result 
holds for strongly non-hermitian matrices, i.e., for Re(z) ~ Im(z) on aver- 
age. In the regime of weak non-hermiticity [14], where the typical magnitude 
of the imaginary parts of the eigenvalues is equal to the mean spacing of the 
real parts, the RMT prediction deviates from Eq. (3). We shall comment 
on this regime below. For uncorrelated eigenvalues in the complex plane, 
the Poisson distribution becomes [11] 



This should not be confused with the Wigner distribution (1). 

Our simulations were done with gauge group SU(3) on a 6 3 x 4 lattice 
using (3 = 6/g 2 = 5.2 in the confinement region and (3 = 5.4 in the decon- 
finement region for Nf = 3 flavors of staggered fermions of mass ma = 0.1. 
Despite major efforts [15] there is currently no feasible solution to the prob- 
lem of a complex weight function in lattice simulations. (In a random matrix 
model, the numerical effort to generate a statistically significant ensemble 
of configurations including the complex Dirac determinant was shown to 
grow exponentially with [i 2 N, where N is the lattice size [16].) Therefore, 
the gauge field configurations were generated at ji = 0, and the chemical po- 
tential was added to the Dirac matrix afterwards. Both in the confinement 
and deconfinement, we sampled 50 independent configurations. 

Typical eigenvalue spectra are shown in Fig. 3 for four different values 
of n (in units of 1/a) at (3 = 5.2. As expected, the size of the real parts of 
the eigenvalues grows with fi, consistent with Ref. [17]. Since the average 
spectral density is not constant, we have to apply the unfolding method 
defined in [12]. 

Our results for P(s) are presented in Fig. 4. There are minor quantita- 
tive but no qualitative differences between confinement and deconfinement 
phase, which is consistent with our findings at \x = (see Fig. 1). As a 
function of /i, we expect to find a transition from Wigner to Ginibre be- 
havior in P(s), as is indeed seen in the figures. For ji = 0.1, the data are 




YN-l 



1 N-l 



(3) 




(4) 



6 



B.A. BERG ET AL. 




Figure 3. Scatter plot of the eigenvalues of the Dirac operator (in units of 1 /a) in the 
complex plane at various values of /i for a typical configuration of full QCD (generated 
at fi = 0) in the confinement region at (3 = 5.2. 

still very close to the Wigner distribution (1) whereas for 0.5 < fi < 0.7 
(fi = 0.7 not shown) we observe nice agreement with the Ginibre distribu- 
tion (3). Values of \x in the crossover region between Wigner and Ginibre 
behavior (0.1 < fi < 0.3) correspond to the regime of weak non-hermiticity 
mentioned above (the plots for \i = 0.3 can be found in Ref. [12]). In this 
regime, the derivation of the spacing distribution is a very difficult prob- 
lem, and the only known analytical result is P(s,zq) for small s, where 
zq is the location in the complex plane (i.e., no unfolding is performed) 
[14]. The small-s behavior of Eqs. (1) and (3) is given by Pw( s ) s 2 an d 
Pq( s ) ^ s 3 ) respectively, and in the regime of weak non-hermiticity we have 
P(s, zq) oc s a (for s«l) with 2 < a < 3 [14]. This smooth crossover from 
oi = 1 to ft = 3 is also observed in our unfolded data. 

For fj, > 0.7 the lattice results for P(s) deviate substantially from the 
Ginibre distribution. The global spectral density of the lattice data for 
fi = 1.0 and 2.2 in Fig. 3 is very different from that of the Ginibre ensem- 
ble. This does not immediately imply that the local spectral fluctuations are 
also different, but it is an indication for qualitative changes. The results for 
fx = 2.2 in Fig. 4 could be interpreted as Poisson behavior, corresponding 
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Figure 4- Nearest-neighbor spacing distribution of the Dirac operator eigenvalues in the 
complex plane for various values of fi in the confinement (left) and deconfinement (right) 
phase. The histograms represent the lattice QCD data. The solid curve is the Ginibrc 
distribution of Eq. (3), the short-dashed curve in the first row the Wigner distribution 
of Eq. (1), and the dotted curve in the last row the Poisson distribution of Eq. (4). 



to uncorrelated eigenvalues. (In the hermitian case at nonzero temperature, 
lattice simulations only show a transition to Poisson behavior for (5 — > oo 
when the physical box size shrinks and the theory becomes free [6].) A 
plausible explanation of the transition to Poisson behavior is provided by 
the following two (related) observations. First, for large \i the terms con- 
taining e M in Eq. (2) dominate the Dirac matrix, giving rise to uncorrelated 
eigenvalues. Second, for fi > 1.0 the fermion density on the 6 3 x 4 lattice 
reaches saturation due to limited box size and the Pauli exclusion principle. 
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Figure 5. Nearest-neighbor spacing distribution P(s) for U(l) gauge theory on an 8 3 x 6 
lattice in the confined phase (left) and in the Coulomb phase (right). The theoretical 
curves are the chUE result, Eq. (1), and the Poisson distribution, Pp(s) = exp( — s). 

3. QED at Nonzero Temperature 

By now it is a well-known fact that the spectrum of the QCD Dirac operator 
is related to universality classes of RMT, i.e., determined by the global 
symmetries of the QCD partition function. We have investigated 4d U(l) 
gauge theory which was not classified yet. At (3 C « 1.01 U(l) gauge theory 
undergoes a phase transition between a confinement phase with mass gap 
and monopole excitations for (3 < (3 C and the Coulomb phase which exhibits 
a massless photon [18] for (3 > (3 C . As for SU(2) and SU(3) gauge groups, we 
expect the confined phase to be described by RMT, whereas free fermions 
are known to yield the Poisson distribution (see Fig. 2). The question arose 
whether the Coulomb phase will be described by RMT or by the Poisson 
distribution [19]. The nearest-neighbor spacing distributions for an 8 3 x 
6 lattice at f3 = 0.9 (confined phase) and at {3 = 1.1 (Coulomb phase), 
averaged over 20 independent configuration, are depicted in Fig. 5. Both 
are well described by the chiral unitary ensemble (chUE) of RMT. 

We have continued the above investigation with a study of the dis- 
tribution of small eigenvalues in the confined phase. The Banks-Casher 
formula [1] relates the eigenvalue density p(\) at A = to the chiral con- 
densate, £ = K^V 7 )! = h m m^o lhnv->oo fl7?(0)/V. The microscopic spectral 
density, p s (z) = limy^oo p (z/VTi) /VS, should be given by the result for 
the chUE of RMT [20]. This function also generates the Leutwyler-Smilga 
sum rules [21]. 

To study the smallest eigenvalues, spectral averaging is not possible, and 
one has to produce large numbers of configurations. Our present results are 
for (3 = 0.9 in the confined phase with 10000 configurations on a 4 4 , 10000 
configuration on a 6 4 , and 2822 configurations on an 8 3 x 6 lattice. The left 
plot in Fig. 6 exhibits the distribution P(A m i n ) of the smallest eigenvalue 
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Figure 6. Distribution P(A m i n ) (left plot) and microscopic spectral density p s (z) (right 
plot) from our 6 4 lattice data of U(l) gauge theory in comparison with the predictions 
of the chUE of RMT (dashed lines). 

A m i n in comparison with the prediction of the (quenched) chUE of RMT 
for topological charge v = 0, 



The agreement is excellent for all lattices. For the chiral condensate we 
obtain S ~ 0.35 by extrapolating the histogram for p(\) to A = and 
using the Banks-Casher relation. Since the average value of A m i n goes like 
V -1 , (A m i n ) decreases with increasing lattice size. In the right plot of Fig. 6 
the same comparison with RMT is done for the microscopic spectral density 
p s (z) up to z = 10, and the agreement is again quite satisfactory. Here, the 
analytical RMT result for the (quenched) chUE and v = is given by [20] 
p s {z) = z [Jo (z) + Jf (z)]/2, where J denotes the Bessel function. 

The quasi-zero modes which are responsible for the chiral condensate 
X i=3 0.35 build up when we cross from the Coulomb into the confined phase. 
For our 8 3 x 6 lattice, Fig. 7 compares on identical scales densities of the 
small eigenvalues at (3 = 0.9 (left plot) and at f3 = 1.1 (right plot), averaged 
over 20 configurations. The quasi-zero modes in the left plot are related 
to the nonzero chiral condensate S > 0, whereas no such quasi-zero modes 
are found in the Coulomb phase. It may be worthwhile to understand the 
physical origin of the U(l) quasi-zero modes in more detail. For 4d SU(2) 
and SU(3) gauge theories a general interpretation is to link them, and 
hence the chiral condensate, to the existence of instantons. As there are 
no instantons in 4d U(l) gauge theory, one needs another explanation, and 
it is interesting to study local correlations of the fermion density with the 
topological charge density and the monopole density [22]. 





10 B.A. BERG ET AL. 




0.1 0.2 ^ 0.1 0.2 ^ 



Figure 7. Density p(A) of small eigenvalues for the 8x6 lattice at f3 = 0.9 (left plot) 
and at /3 = 1.1 (right plot). A nonzero chiral condensate is supported in the confinement 
phase of U(l) gauge theory. 



Another interesting question concerns the energy scale E c up to which 
RMT describes the small Dirac eigenvalues in the phase where p(0) > 0. 
In disordered mesoscopic systems, a similar scale is called the Thouless 
energy. The theoretical prediction for QCD is E c ~ f%/T,L^ [23] with the 
pion decay constant f n , where we have assumed that the spatial extent L s 
of the lattice is not smaller than the temporal extent Lt- In units of the 
mean level spacing A = ir/VT, at the origin, this becomes 

u c = ^ ~ -flL s L t . (6) 

A convenient quantity from which u c can be extracted is the disconnected 
scalar susceptibility, 

The corresponding RMT result for the quenched chUE with v = reads 

[24] Xrmt = u2 i K H u ) ~ K o( u Wo( u ) ~ ^(u)], where u = mV ^> and 1 and 
K are modified Bessel functions. In Fig. 8 we have plotted the ratio [25] 



ratio 



.disc disc \ /^ / disc /o\ 

l^Xiatt ~ Xrmt) /armt \°) 

versus u and uj (L s L t ), respectively, for the U(l) data computed at (5 = 0.9. 
This ratio should deviate from zero above the Thouless scale. The expected 
scaling of the Thouless energy with L s Lt is confirmed. 



4. Conclusions 

We have searched for a transition in the nearest-neighbor spacing distribu- 
tion P(s) from Wigner to Poisson behavior across the deconfinement phase 
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Figure 8. The ratio of Eq. (8) for U(l) gauge theory plotted versus u and u/(L s Lt), 
respectively (error bars not shown). In the right plot, the data for different L a and L t fall 
on the same curve, confirming the expected scaling of the Thouless energy according to 
Eq. (6) . The deviations of the ratio from zero for very small values of u are well-understood 
artifacts of the finite lattice size and finite statistics [25]. 

transition of pure gluonic and of full QCD. We observed no signature of a 
transition, neither for pure SU(3) nor for full QCD. The data agree with the 
RMT result in both phases, except for extremely large values of (5 where 
the eigenvalues are known analytically. Our analysis of full QCD shows that 
quark masses have no influence on the nearest-neighbor spacing distribu- 
tion. One explanation of our results is that temporal monopole currents 
survive the deconfinement phase transition leading to confinement of spa- 
tial Wilson loops. Thus, even in the deconfinement phase, the gauge fields 
retain a considerable degree of randomness. 

A general unfolding procedure for the spectra of non-hermitian oper- 
ators was applied to the QCD lattice Dirac operator at nonzero chemical 
potential. Agreement of the nearest-neighbor spacing distribution with pre- 
dictions of the Ginibre ensemble of non-hermitian RMT was found between 
/i = 0.5 and ji = 0.7 in both confinement and deconfinement phase. The de- 
viations from Ginibre behavior for smaller values of \i as well as the changes 
for larger values of \i toward a Poisson shape are understood mathemati- 
cally. The physical interpretation requires a better understanding of QCD 
at nonzero density. An interesting observation is that the results for P(s) 
in the non-hermitian case are rather sensitive to [i whereas they are very 
stable under variations of T in the hermitian case. 

The nearest-neighbor spacing distribution of 4d U(l) quenched lattice 
gauge theory is described by the chUE in both the confinement and the 
Coulomb phase. In the confinement phase we also find that the -P(A m i n ) 
distribution and the microscopic spectral density p s (z) are described by 
the chUE. The Thouless energy scales with the lattice size as expected. 



12 



B.A. BERG ET AL. 



5. Acknowledgments 

This work was supported in part by FWF projects P10468-PHY and P11456- 
PHY, by DFG grants We 655/11-2 and We 655/15-1, by DOE contracts DE- 
FG02-97ER41022, DE-FG05-85ER2500, and DE-FG02-91ER40608, and by 
the RIKEN BNL Research Center. We thank T.S. Biro, E.-M. Ilgenfritz, N. 
Kaiser, M.I. Polikarpov, K. Rabitsch, and J.J.M. Verbaarschot for helpful 
discussions. 



References 

1. T. Banks and A. Casher, Nucl. Phys. B 169 (1980) 103. 

2. M.A. Halasz and J.J.M. Verbaarschot, Phys. Rev. Lett. 74 (1995) 3920; M.A. Halasz, 
T. Kalkreuter, and J.J.M. Verbaarschot, Nucl. Phys. B (Proc. Suppl.) 53 (1997) 266. 

3. O. Bohigas, M.-J. Giannoni, and C. Schmit, Phys. Rev. Lett. 52 (1984) 1. 

4. J.J.M. Verbaarschot, Phys. Rev. Lett. 72 (1994) 2531. 

5. D. Fox and P.B. Kahn, Phys. Rev. 134 (1964) B1151; T. Nagao and M. Wadati, J. 
Phys. Soc. Jpn. 60 (1991) 3298; 61 (1992) 78; 61 (1992) 1910. 

6. R. Pullirsch et al., Phys. Lett. B 427 (1998) 119. 

7. M.L. Mehta, Random Matrices, 2nd ed. (Academic Press, San Diego, 1991). 

8. M.A. Stephanov, Phys. Rev. Lett. 76 (1996) 4472. 

9. P. Hasenfratz and F. Karsch, Phys. Lett. B 125 (1983) 308; I.M. Barbour, Nucl. 
Phys. B (Proc. Suppl.) 26 (1992) 22. 

10. M.A. Halasz, J.C. Osborn, and J.J.M. Verbaarschot, Phys. Rev. D 56 (1997) 7059. 

11. R. Grobe, F. Haake, and H.-J. Sommers, Phys. Rev. Lett. 61 (1988) 1899. 

12. H. Markum, R. Pullirsch, and T. Wettig, Phys. Rev. Lett. 83 (1999) 484. 

13. J. Ginibre, J. Math. Phys. 6 (1965) 440. 

14. Y.V. Fyodorov and H.-J. Sommers, JETP Lett. 63 (1996) 1026; J. Math. Phys. 38 
(1997) 1918; Y.V. Fyodorov, B.A. Khoruzhenko, and H.-J. Sommers, Phys. Lett. A 
226 (1997) 46; Phys. Rev. Lett. 79 (1997) 557. 

15. I.M. Barbour et al., Phys. Rev. D 56 (1997) 7063. 

16. M.A. Halasz, A.D. Jackson, and J.J.M. Verbaarschot, Phys. Rev. D 56 (1997) 5140. 

17. I. Barbour et al, Nucl. Phys. B 275 (1986) 296. 

18. B.A. Berg and C. Panagiotakopoulos, Phys. Rev. Lett. 52 (1984) 94. 

19. B.A. Berg, H. Markum, and R. Pullirsch, Phys. Rev. D 59 (1999) 097504. 

20. E.V. Shuryak and J.J.M. Verbaarschot, Nucl. Phys. A 560 (1992) 306; J.J.M. Ver- 
baarschot and I. Zahed, Phys. Rev. Lett. 70 (1993) 3852; J.C. Osborn, D. Toublan, 
and J.J.M. Verbaarschot, Nucl. Phys. B 540 (1999) 317; P.H. Damgaard et al., Nucl. 
Phys. B 547 (1999) 305. 

21. H. Leutwyler and A.V. Smilga, Phys. Rev. D 46 (1992) 5607. 

22. W. Sakuler, S. Thurner, and H. Markum, Phys. Lett. B 464 (1999) 272. 

23. J.C. Osborn and J.J.M. Verbaarschot, Nucl. Phys. B 525 (1998) 738; Phys. Rev. 
Lett. 81 (1998) 268; R.A. Janik et al., Phys. Rev. Lett. 81 (1998) 264. 

24. M. Gockeler et al., Phys. Rev. D 59 (1999) 094503. 

25. M.E. Berbenni-Bitsch et al., Phys. Lett. B 438 (1998) 14. 



